NOTES FROM SF:

FOR STEAMBOAT, COMPARE THE DETRENDED CLIAMTE NORMAL CURVES AGAINST EACH OTHER

Derry & fassnacht 2010 Compare climatology of snotel statuions= they don’t have the same climatol9ogy based on the snow data

Figure 5b - https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2009WR007835

Disclaimer

This work is very preliminary as I get back into the coding swing of things. Data wrangling and figure generation will be done via R, but the rest of the project will be done using good ol’ microsoft products. This is just an entry point into data crunching and should by no means be considered a final product.

Also, I’m not great at this but whatever. I could automate this, but I’ll figure that out shortly!

Need to update figures & analyze sites for seasonal trends

Methodology

SNOTEL data was provided by the NRCS. Data was cleaned by removing outliers that are likely implausible; any year with more than 15 observations missing was removed. Temperatures were adjusted using the Morrisey method for stations identified by Ma et al (2019) due to SNOTEL temperature sensor changes, with the adjustment applied to pre-sensor change data. Daily mean observations were detrended to determine whether values were increasing or decreasing from the entire time series trend. Daily mean temperatures were first averaged by water year, with all water year means then averaged by day of water year. The mean temperature by day for the period of record was averaged. To find the standard deviation, the daily mean temperatures by water year was subtracted from the averaged mean temperature by day for the period of record. All water year means averaged by day of water year were subtracted from the temperature mean. The resulting values were then added together to find the “residual” of the daily mean temperatures by water year. The standard deviation was then computed from those residuals, with trends analyzed by Mann‐Kendall significance test and Theil‐Sen’s rate of change. Significant trends are identified with p-values of less than 0.10.

Morrisey Method

The Morrisey Method is taken from Ma, Fassnacht and Kampf..

In R script: T(adjusted) = 5.3x10(-7)xT(old)4+3.72x10(-5)xT(old)3-2.16x10(-3)xT(old)2-7.32x10^(-2)xT(old)+1.37

In the Ma et al. spreadsheet, H1 is Morrisey, H2 is Oiler

Yampa Area SNOTEL sites:

Burro Mountain 378 NSCE- Original, Bias- Morrisey 7/11/2002 (NSCE- 0.87 vs 0.84, Bias- -0.03 vs. 0.11)

Columbine 408 Morrisey 7/22/2005

Columbine Pass 409 Morrisey 6/23/2005

Crosho 426 Morrisey 7/21/2005

Dry Lake 457 Oyler -> Morrisey 7/30/2003

Elk River 467 Oyler -> Morrisey 8/7/2006

Lynx Pass 607 Morrisey 5/22/2006

Rabbit Ears 709 Morrisey 8/7/2006 (Does not go.)

Ripple Creek 717 Oyler -> Morrisey 8/7/2006

Tower 825 Original 8/18/2004

Trapper Lake 827 Original 12/13/2004

Data from thesis research:


SNOTEL_yampa_area <- snotel_download(site_id = c(378, 408, 409, 426, 457, 467, 607, 709, 717, 825, 827), path = tempdir('../data'), internal = TRUE)

write.csv(SNOTEL_yampa_area,"C:/Users/13074/Documents/ESS580/thesis_project/Yampa/data_raw/snotel_yampa.csv", row.names = FALSE) #write in the raw data
SNOTEL_yampa_area <- read.csv("C:/Users/13074/Documents/ESS580/thesis_project/Yampa/data_raw/snotel_yampa.csv", header = TRUE)

Burro Mountain 378

NSCE- Original, Bias- Morrisey 7/11/2002 **Original is marked as having the smallest bias, but the sheet says original is smaller.)

snotel_378 <- SNOTEL_yampa_area %>% 
  filter(site_id == "378")
#str(snotel_378) # check the date, usually a character.  

snotel_378$Date <- as.Date(snotel_378$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.

#THIS WILL CHANGE FOR EACH STATION
snotel_378_clean <- snotel_378 %>% # filter for the timeframe
  filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
  #filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers   
  addWaterYear() %>% 
  mutate(daymonth = format(as.Date(Date), "%d-%m")) %>% 
  na.omit()

#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))

snotel_378_clean <- snotel_378_clean %>% 
  group_by(waterYear)%>% 
  mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers

ggplot(snotel_378_clean, aes(x = Date, y = temperature_mean)) +
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

snotel_378_clean <- snotel_378_clean %>% 
  mutate(temp_diff = abs(temperature_min - temperature_max)) %>% 
  filter(temperature_mean > -40)# %>% 
  #filter(temperature_mean < 25)
ggplot(snotel_378_clean, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”

# filtering for temperature anomalies
#snotel_378_cull_count <- snotel_378_clean %>% 
#  filter(temperature_min > -40) %>% 
#  count(waterYear)

#snotel_378_cull_count

# filtering for too few observations in a year
snotel_378_cull_count_days <- snotel_378_clean %>% 
  group_by(waterYear) %>% 
  count(waterYear) %>% 
  filter(n < 350)

snotel_378_cull_count_days
## # A tibble: 10 x 2
## # Groups:   waterYear [10]
##    waterYear     n
##        <dbl> <int>
##  1      1985     1
##  2      1988   298
##  3      1989   341
##  4      1993   281
##  5      1994   261
##  6      1995   344
##  7      1996   320
##  8      2002   312
##  9      2016   310
## 10      2022   338
snotel_378_clean_culled <- snotel_378_clean %>% 
  filter(waterYear != "1985" & waterYear != "1986" & waterYear != "1988" & waterYear != "1989" & waterYear != "1993" & waterYear != "1994" & waterYear != "1995" & waterYear != "1996" & waterYear != "2002" & waterYear != "2016" & waterYear != "2022")# & waterYear != "2017") #%>% 
  #filter(temperature_mean > -49)

ggplot(snotel_378_clean_culled, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

Culling WY 1986 as well.

ggplot(snotel_378_clean_culled, aes(x = Date, y = temperature_mean)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

temp_378_xts <- xts(snotel_378_clean_culled$temperature_mean, order.by = snotel_378_clean_culled$Date)

dygraph(temp_378_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 
#snotel_378_clean_culled <- snotel_378_clean_culled %>% 
#  filter(temperature_mean < 19.3)

temp_378_xts <- xts(snotel_378_clean_culled$temperature_mean, order.by = snotel_378_clean_culled$Date)

dygraph(temp_378_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 

Burro Mountain 378

NSCE- Original, Bias- Morrisey 7/11/2002

snotel_378_adjusted <- snotel_378_clean_culled %>%
  mutate(temp_ad = if_else(Date < "2002-07-11", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))

Non-corrected

378 Detrended

#using the clean culled df:

#average water year temperature

yearly_wy_aver_378 <- snotel_378_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(temperature_mean))

#Average temperature by day for all water years:

daily_wy_aver_378 <- yearly_wy_aver_378 %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

daily_wy_aver_378 <- daily_wy_aver_378 %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(daily_wy_aver_378$aver_day_temp))

# try to show all years as means. 
daily_wy_aver2_378 <-daily_wy_aver_378 %>% 
  group_by(waterDay) %>%
  mutate(date_temp = mean(temperature_mean))
  

daily_wy_aver2_378$date_temp <- signif(daily_wy_aver2_378$date_temp,3) #reduce the sig figs


ggplot(daily_wy_aver2_378, aes(x = waterDay, y = date_temp))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

378 SD

standard_dev_378 <- daily_wy_aver_378 %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

standard_dev_all_378 <- standard_dev_378 %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

standard_dev_all_378 <- standard_dev_all_378 %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

standard_dev_all_378 %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1987 8.562568
1990 9.114355
1991 9.254156
1992 8.275786
1997 8.898932
1998 8.983050
1999 8.073327
2000 8.686540
2001 9.301471
2003 8.509538
2004 8.017044
2005 7.723002
2006 8.547829
2007 8.796189
2008 8.749038
2009 7.771973
2010 8.617066
2011 8.282002
2012 8.445614
2013 8.991985
2014 8.166025
2015 7.657629
2017 8.034586
2018 8.020621
2019 8.518154
2020 8.848254
2021 8.718771
ggplot(standard_dev_all_378, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Non-corrected standard deviation of SNOTEL 378 average temperatures for water years 2005-2021

MK & SS for 378 (non-corrected)

sd_mk_378 <- mk.test(standard_dev_all_378$sd_2)
print(sd_mk_378)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_378$sd_2
## z = -1.5427, n = 27, p-value = 0.1229
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  -75.0000000 2301.0000000   -0.2136752
sd_sens_378 <- sens.slope(standard_dev_all_378$sd_2)
print(sd_sens_378)
## 
##  Sen's slope
## 
## data:  standard_dev_all_378$sd_2
## z = -1.5427, n = 27, p-value = 0.1229
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.047456119  0.006621254
## sample estimates:
## Sen's slope 
##  -0.0230883

Corrected

378 Detrended (corrected)

#using the clean culled df:
#average water year temperature
yearly_wy_aver_378_ad <- snotel_378_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_378_ad <- yearly_wy_aver_378_ad %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_378_ad <- daily_wy_aver_378_ad %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp_ad = mean(daily_wy_aver_378_ad$aver_day_temp_ad))
# try to show all years as means. 
daily_wy_aver2_378_ad <-daily_wy_aver_378_ad %>% 
  group_by(waterDay) %>%
  mutate(date_temp_ad = mean(temp_ad))
  
daily_wy_aver2_378_ad$date_temp_ad <- signif(daily_wy_aver2_378_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_378_ad, aes(x = waterDay, y = date_temp_ad))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

378 SS (corrected)

standard_dev_378_ad <- daily_wy_aver_378_ad %>% 
  group_by(waterYear) %>% 
  #filter(waterYear >= 1987 & waterYear <= 2021) %>% 
  mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>% 
  mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_378_ad <- standard_dev_378_ad %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())
standard_dev_all_378_ad <- standard_dev_all_378_ad %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)
standard_dev_all_378_ad %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1987 7.940547
1990 8.436056
1991 8.633602
1992 7.670874
1997 8.260208
1998 8.299486
1999 7.471734
2000 7.993287
2001 8.584377
2003 8.508522
2004 8.016141
2005 7.722598
2006 8.547499
2007 8.795362
2008 8.748133
2009 7.771479
2010 8.616338
2011 8.281602
2012 8.444783
2013 8.991116
2014 8.165236
2015 7.657576
2017 8.034186
2018 8.019986
2019 8.517552
2020 8.847387
2021 8.717958
ggplot(standard_dev_all_378_ad, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Morrisey-corrected standard deviation of SNOTEL 378 average temperatures for water years 1986-2021

MK & SS 378 (corrected)

sd_mk_378_ad <- mk.test(standard_dev_all_378_ad$sd_2)
print(sd_mk_378_ad)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_378_ad$sd_2
## z = 1.3342, n = 27, p-value = 0.1821
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##   65.0000000 2301.0000000    0.1851852
sd_sens_378_ad <- sens.slope(standard_dev_all_378_ad$sd_2)
print(sd_sens_378_ad)
## 
##  Sen's slope
## 
## data:  standard_dev_all_378_ad$sd_2
## z = 1.3342, n = 27, p-value = 0.1821
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.009653315  0.035911433
## sample estimates:
## Sen's slope 
##  0.01222034

Columbine 408

Morrisey 7/22/2005

snotel_408 <- SNOTEL_yampa_area %>% 
  filter(site_id == "408")
#str(snotel_408) # check the date, usually a character.  

snotel_408$Date <- as.Date(snotel_408$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.

#THIS WILL CHANGE FOR EACH STATION
snotel_408_clean <- snotel_408 %>% # filter for the timeframe
  filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
  #filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers   
  addWaterYear() %>% 
  mutate(daymonth = format(as.Date(Date), "%d-%m")) %>% 
  na.omit()

#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))

snotel_408_clean <- snotel_408_clean %>% 
  group_by(waterYear)%>% 
  mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers

ggplot(snotel_408_clean, aes(x = Date, y = temperature_mean)) +
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

snotel_408_clean <- snotel_408_clean %>% 
  mutate(temp_diff = abs(temperature_min - temperature_max)) %>% 
  filter(temperature_mean > -40)# %>% 
  #filter(temperature_mean < 25)
ggplot(snotel_408_clean, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”

# filtering for temperature anomalies
#snotel_408_cull_count <- snotel_408_clean %>% 
#  filter(temperature_min > -40) %>% 
#  count(waterYear)

#snotel_408_cull_count

# filtering for too few observations in a year
snotel_408_cull_count_days <- snotel_408_clean %>% 
  group_by(waterYear) %>% 
  count(waterYear) %>% 
  filter(n < 350)

snotel_408_cull_count_days
## # A tibble: 8 x 2
## # Groups:   waterYear [8]
##   waterYear     n
##       <dbl> <int>
## 1      1981   343
## 2      1983   317
## 3      1984   234
## 4      1985   195
## 5      1986   336
## 6      1987   298
## 7      1994   272
## 8      2002   347
snotel_408_clean_culled <- snotel_408_clean %>% 
  filter(waterYear > "1987" & waterYear != "1993" & waterYear != "1994" & waterYear != "2002")# & waterYear != "1985" & waterYear != "1986" & waterYear != "1987" & waterYear != "1994" & waterYear != "2002")# & waterYear != "2002" & waterYear != "2016" & waterYear != "2022")# & waterYear != "2017") #%>% 
  #filter(temperature_mean > -49)

ggplot(snotel_408_clean_culled, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

Culling WY 1986 & 1993 as well.

ggplot(snotel_408_clean_culled, aes(x = Date, y = temperature_mean)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

temp_408_xts <- xts(snotel_408_clean_culled$temperature_mean, order.by = snotel_408_clean_culled$Date)

dygraph(temp_408_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 
snotel_408_clean_culled <- snotel_408_clean_culled %>% 
  filter(temperature_mean > -30)

temp_408_xts <- xts(snotel_408_clean_culled$temperature_mean, order.by = snotel_408_clean_culled$Date)

dygraph(temp_408_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 

Columbine 408 Morrisey 7/22/2005

snotel_408_adjusted <- snotel_408_clean_culled %>%
  mutate(temp_ad = if_else(Date < "2005-07-22", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))

Non-corrected

408 Detrended

#using the clean culled df:

#average water year temperature

yearly_wy_aver_408 <- snotel_408_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(temperature_mean))

#Average temperature by day for all water years:

daily_wy_aver_408 <- yearly_wy_aver_408 %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

daily_wy_aver_408 <- daily_wy_aver_408 %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(daily_wy_aver_408$aver_day_temp))

# try to show all years as means. 
daily_wy_aver2_408 <-daily_wy_aver_408 %>% 
  group_by(waterDay) %>%
  mutate(date_temp = mean(temperature_mean))
  

daily_wy_aver2_408$date_temp <- signif(daily_wy_aver2_408$date_temp,3) #reduce the sig figs


ggplot(daily_wy_aver2_408, aes(x = waterDay, y = date_temp))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

408 SD

standard_dev_408 <- daily_wy_aver_408 %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

standard_dev_all_408 <- standard_dev_408 %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

standard_dev_all_408 <- standard_dev_all_408 %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

standard_dev_all_408 %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1988 9.911311
1989 9.479817
1990 8.845588
1991 9.141713
1992 8.235717
1995 8.387811
1996 8.807923
1997 8.867338
1998 9.074387
1999 8.232888
2000 8.585091
2001 9.380546
2003 8.847298
2004 8.495495
2005 8.135796
2006 8.734304
2007 8.839901
2008 8.888992
2009 7.926154
2010 8.738675
2011 8.455577
2012 8.702054
2013 9.110505
2014 8.268271
2015 7.870630
2016 8.492743
2017 8.238399
2018 8.140089
2019 8.644103
2020 9.027143
2021 8.902353
2022 8.706784
ggplot(standard_dev_all_408, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Non-corrected standard deviation of SNOTEL 408 average temperatures for water years 2005-2021

MK & SS for 408 (non-corrected)

sd_mk_408 <- mk.test(standard_dev_all_408$sd_2)
print(sd_mk_408)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_408$sd_2
## z = -1.7676, n = 32, p-value = 0.07713
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## -110.0000000 3802.6666667   -0.2217742
sd_sens_408 <- sens.slope(standard_dev_all_408$sd_2)
print(sd_sens_408)
## 
##  Sen's slope
## 
## data:  standard_dev_all_408$sd_2
## z = -1.7676, n = 32, p-value = 0.07713
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.033742130  0.001027723
## sample estimates:
## Sen's slope 
## -0.01615257

Corrected

408 Detrended (corrected)

#using the clean culled df:
#average water year temperature
yearly_wy_aver_408_ad <- snotel_408_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_408_ad <- yearly_wy_aver_408_ad %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_408_ad <- daily_wy_aver_408_ad %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp_ad = mean(daily_wy_aver_408_ad$aver_day_temp_ad))
# try to show all years as means. 
daily_wy_aver2_408_ad <-daily_wy_aver_408_ad %>% 
  group_by(waterDay) %>%
  mutate(date_temp_ad = mean(temp_ad))
  
daily_wy_aver2_408_ad$date_temp_ad <- signif(daily_wy_aver2_408_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_408_ad, aes(x = waterDay, y = date_temp_ad))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

408 SS (corrected)

standard_dev_408_ad <- daily_wy_aver_408_ad %>% 
  group_by(waterYear) %>% 
  #filter(waterYear >= 1987 & waterYear <= 2021) %>% 
  mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>% 
  mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_408_ad <- standard_dev_408_ad %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())
standard_dev_all_408_ad <- standard_dev_all_408_ad %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)
standard_dev_all_408_ad %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1988 9.255004
1989 8.873391
1990 8.216764
1991 8.545580
1992 7.664569
1995 7.780528
1996 8.186934
1997 8.260524
1998 8.410419
1999 7.651225
2000 7.928801
2001 8.691699
2003 8.174582
2004 7.895518
2005 7.482674
2006 8.735962
2007 8.841502
2008 8.890875
2009 7.928502
2010 8.741101
2011 8.457260
2012 8.703940
2013 9.112301
2014 8.270874
2015 7.872234
2016 8.495258
2017 8.239886
2018 8.142322
2019 8.646452
2020 9.028953
2021 8.904182
2022 8.708980
ggplot(standard_dev_all_408_ad, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Morrisey-corrected standard deviation of SNOTEL 408 average temperatures for water years 1986-2021

MK & SS 408 (corrected)

sd_mk_408_ad <- mk.test(standard_dev_all_408_ad$sd_2)
print(sd_mk_408_ad)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_408_ad$sd_2
## z = 1.2487, n = 32, p-value = 0.2118
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##   78.0000000 3802.6666667    0.1572581
sd_sens_408_ad <- sens.slope(standard_dev_all_408_ad$sd_2)
print(sd_sens_408_ad)
## 
##  Sen's slope
## 
## data:  standard_dev_all_408_ad$sd_2
## z = 1.2487, n = 32, p-value = 0.2118
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.007349474  0.033195852
## sample estimates:
## Sen's slope 
##  0.01316118

Columbine Pass 409

Morrisey 6/23/2005

snotel_409 <- SNOTEL_yampa_area %>% 
  filter(site_id == "409")
#str(snotel_409) # check the date, usually a character.  

snotel_409$Date <- as.Date(snotel_409$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.

#THIS WILL CHANGE FOR EACH STATION
snotel_409_clean <- snotel_409 %>% # filter for the timeframe
  filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
  #filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers   
  addWaterYear() %>% 
  mutate(daymonth = format(as.Date(Date), "%d-%m")) %>% 
  na.omit()

#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))

snotel_409_clean <- snotel_409_clean %>% 
  group_by(waterYear)%>% 
  mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers

ggplot(snotel_409_clean, aes(x = Date, y = temperature_mean)) +
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

snotel_409_clean <- snotel_409_clean %>% 
  mutate(temp_diff = abs(temperature_min - temperature_max)) %>% 
  filter(temperature_mean > -50) %>% 
  filter(temperature_mean < 40)
ggplot(snotel_409_clean, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”

# filtering for temperature anomalies
#snotel_409_cull_count <- snotel_409_clean %>% 
#  filter(temperature_min > -40) %>% 
#  count(waterYear)

#snotel_409_cull_count

# filtering for too few observations in a year
snotel_409_cull_count_days <- snotel_409_clean %>% 
  group_by(waterYear) %>% 
  count(waterYear) %>% 
  filter(n < 350)

snotel_409_cull_count_days
## # A tibble: 5 x 2
## # Groups:   waterYear [5]
##   waterYear     n
##       <dbl> <int>
## 1      1994   331
## 2      1995   190
## 3      2005   329
## 4      2019   348
## 5      2022   346
snotel_409_clean_culled <- snotel_409_clean %>% 
  filter(waterYear != "1994" & waterYear != "1995" & waterYear != "2005" & waterYear != "2019" & waterYear != "2022")# & waterYear != "1986" & waterYear != "1987" & waterYear != "1994" & waterYear != "2002")# & waterYear != "2002" & waterYear != "2016" & waterYear != "2022")# & waterYear != "2017") #%>% 
  #filter(temperature_mean > -49)

ggplot(snotel_409_clean_culled, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

ggplot(snotel_409_clean_culled, aes(x = Date, y = temperature_mean)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

temp_409_xts <- xts(snotel_409_clean_culled$temperature_mean, order.by = snotel_409_clean_culled$Date)

dygraph(temp_409_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 
#snotel_409_clean_culled <- snotel_409_clean_culled %>% 
#  filter(temperature_mean > -30)

#temp_409_xts <- xts(snotel_409_clean_culled$temperature_mean, order.by = snotel_409_clean_culled$Date)

#dygraph(temp_409_xts) %>%
#  dyAxis("y", label = "Daily mean temperature (°C)") 

Columbine Pass 409

Morrisey 6/23/2005

snotel_409_adjusted <- snotel_409_clean_culled %>%
  mutate(temp_ad = if_else(Date < "2005-06-23", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))

Non-corrected

409 Detrended

#using the clean culled df:

#average water year temperature

yearly_wy_aver_409 <- snotel_409_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(temperature_mean))

#Average temperature by day for all water years:

daily_wy_aver_409 <- yearly_wy_aver_409 %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

daily_wy_aver_409 <- daily_wy_aver_409 %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(daily_wy_aver_409$aver_day_temp))

# try to show all years as means. 
daily_wy_aver2_409 <-daily_wy_aver_409 %>% 
  group_by(waterDay) %>%
  mutate(date_temp = mean(temperature_mean))
  

daily_wy_aver2_409$date_temp <- signif(daily_wy_aver2_409$date_temp,3) #reduce the sig figs


ggplot(daily_wy_aver2_409, aes(x = waterDay, y = date_temp))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

409 SD

standard_dev_409 <- daily_wy_aver_409 %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

standard_dev_all_409 <- standard_dev_409 %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

standard_dev_all_409 <- standard_dev_all_409 %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

standard_dev_all_409 %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1987 8.326403
1988 9.167257
1989 8.983106
1990 8.555010
1991 7.839925
1992 7.987351
1993 8.471182
1996 8.508682
1997 8.478813
1998 8.929006
1999 7.623863
2000 8.631812
2001 9.047367
2002 9.516272
2003 8.791318
2004 8.708583
2006 8.125432
2007 8.536611
2008 8.816576
2009 7.846243
2010 8.511418
2011 8.496243
2012 8.285134
2013 8.792304
2014 7.958185
2015 7.393468
2016 8.187794
2017 8.066498
2018 7.805356
2020 8.856536
2021 8.674169
ggplot(standard_dev_all_409, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Non-corrected standard deviation of SNOTEL 409 average temperatures for water years 2005-2021

MK & SS for 409 (non-corrected)

sd_mk_409 <- mk.test(standard_dev_all_409$sd_2)
print(sd_mk_409)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_409$sd_2
## z = -1.1218, n = 31, p-value = 0.262
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##           S        varS         tau 
##  -67.000000 3461.666667   -0.144086
sd_sens_409 <- sens.slope(standard_dev_all_409$sd_2)
print(sd_sens_409)
## 
##  Sen's slope
## 
## data:  standard_dev_all_409$sd_2
## z = -1.1218, n = 31, p-value = 0.262
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.034898600  0.009250748
## sample estimates:
## Sen's slope 
## -0.01256069

Corrected

409 Detrended (corrected)

#using the clean culled df:
#average water year temperature
yearly_wy_aver_409_ad <- snotel_409_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_409_ad <- yearly_wy_aver_409_ad %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_409_ad <- daily_wy_aver_409_ad %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp_ad = mean(daily_wy_aver_409_ad$aver_day_temp_ad))
# try to show all years as means. 
daily_wy_aver2_409_ad <-daily_wy_aver_409_ad %>% 
  group_by(waterDay) %>%
  mutate(date_temp_ad = mean(temp_ad))
  
daily_wy_aver2_409_ad$date_temp_ad <- signif(daily_wy_aver2_409_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_409_ad, aes(x = waterDay, y = date_temp_ad))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

409 SS (corrected)

standard_dev_409_ad <- daily_wy_aver_409_ad %>% 
  group_by(waterYear) %>% 
  #filter(waterYear >= 1987 & waterYear <= 2021) %>% 
  mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>% 
  mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_409_ad <- standard_dev_409_ad %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())
standard_dev_all_409_ad <- standard_dev_all_409_ad %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)
standard_dev_all_409_ad %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1987 7.660579
1988 8.456632
1989 8.306699
1990 7.862004
1991 7.209890
1992 7.340956
1993 7.806529
1996 7.809598
1997 7.810012
1998 8.200421
1999 7.004624
2000 7.910753
2001 8.302261
2002 8.749019
2003 8.046898
2004 8.036896
2006 8.125530
2007 8.536428
2008 8.816453
2009 7.846083
2010 8.511060
2011 8.495994
2012 8.285156
2013 8.792015
2014 7.958138
2015 7.393110
2016 8.187413
2017 8.066215
2018 7.805282
2020 8.856204
2021 8.673923
ggplot(standard_dev_all_409_ad, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Morrisey-corrected standard deviation of SNOTEL 409 average temperatures for water years 1986-2021

MK & SS 409 (corrected)

sd_mk_409_ad <- mk.test(standard_dev_all_409_ad$sd_2)
print(sd_mk_409_ad)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_409_ad$sd_2
## z = 2.1415, n = 31, p-value = 0.03223
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  127.0000000 3461.6666667    0.2731183
sd_sens_409_ad <- sens.slope(standard_dev_all_409_ad$sd_2)
print(sd_sens_409_ad)
## 
##  Sen's slope
## 
## data:  standard_dev_all_409_ad$sd_2
## z = 2.1415, n = 31, p-value = 0.03223
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  0.002607708 0.039439754
## sample estimates:
## Sen's slope 
##  0.02152587

Crosho 426

Morrisey 7/21/2005

snotel_426 <- SNOTEL_yampa_area %>% 
  filter(site_id == "426")
#str(snotel_426) # check the date, usually a character.  

snotel_426$Date <- as.Date(snotel_426$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.

#THIS WILL CHANGE FOR EACH STATION
snotel_426_clean <- snotel_426 %>% # filter for the timeframe
  filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
  #filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers   
  addWaterYear() %>% 
  mutate(daymonth = format(as.Date(Date), "%d-%m")) %>% 
  na.omit()

#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))

snotel_426_clean <- snotel_426_clean %>% 
  group_by(waterYear)%>% 
  mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers

ggplot(snotel_426_clean, aes(x = Date, y = temperature_mean)) +
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

snotel_426_clean <- snotel_426_clean %>% 
  mutate(temp_diff = abs(temperature_min - temperature_max)) %>% 
  filter(temperature_mean > -50) %>% 
  filter(temperature_mean < 40)
ggplot(snotel_426_clean, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”

# filtering for temperature anomalies
#snotel_426_cull_count <- snotel_426_clean %>% 
#  filter(temperature_min > -40) %>% 
#  count(waterYear)

#snotel_426_cull_count

# filtering for too few observations in a year
snotel_426_cull_count_days <- snotel_426_clean %>% 
  group_by(waterYear) %>% 
  count(waterYear) %>% 
  filter(n < 350)

snotel_426_cull_count_days
## # A tibble: 2 x 2
## # Groups:   waterYear [2]
##   waterYear     n
##       <dbl> <int>
## 1      2009   340
## 2      2021   349
snotel_426_clean_culled <- snotel_426_clean %>% 
  filter(waterYear != "2009" & waterYear != "2021")# & waterYear != "2005" & waterYear != "2019" & waterYear != "2022")# & waterYear != "1986" & waterYear != "1987" & waterYear != "1994" & waterYear != "2002")# & waterYear != "2002" & waterYear != "2016" & waterYear != "2022")# & waterYear != "2017") #%>% 
  #filter(temperature_mean > -49)

ggplot(snotel_426_clean_culled, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

ggplot(snotel_426_clean_culled, aes(x = Date, y = temperature_mean)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

temp_426_xts <- xts(snotel_426_clean_culled$temperature_mean, order.by = snotel_426_clean_culled$Date)

dygraph(temp_426_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 
#snotel_426_clean_culled <- snotel_426_clean_culled %>% 
#  filter(temperature_mean > -30)

#temp_426_xts <- xts(snotel_426_clean_culled$temperature_mean, order.by = snotel_426_clean_culled$Date)

#dygraph(temp_426_xts) %>%
#  dyAxis("y", label = "Daily mean temperature (°C)") 

Crosho 426 Morrisey 7/21/2005

snotel_426_adjusted <- snotel_426_clean_culled %>%
  mutate(temp_ad = if_else(Date < "2005-07-21", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))

Non-corrected

426 Detrended

#using the clean culled df:

#average water year temperature

yearly_wy_aver_426 <- snotel_426_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(temperature_mean))

#Average temperature by day for all water years:

daily_wy_aver_426 <- yearly_wy_aver_426 %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

daily_wy_aver_426 <- daily_wy_aver_426 %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(daily_wy_aver_426$aver_day_temp))

# try to show all years as means. 
daily_wy_aver2_426 <-daily_wy_aver_426 %>% 
  group_by(waterDay) %>%
  mutate(date_temp = mean(temperature_mean))
  

daily_wy_aver2_426$date_temp <- signif(daily_wy_aver2_426$date_temp,3) #reduce the sig figs


ggplot(daily_wy_aver2_426, aes(x = waterDay, y = date_temp))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

426 SD

standard_dev_426 <- daily_wy_aver_426 %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

standard_dev_all_426 <- standard_dev_426 %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

standard_dev_all_426 <- standard_dev_all_426 %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

standard_dev_all_426 %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1987 8.642717
1988 8.008775
1989 9.449036
1990 8.683053
1991 9.170577
1992 8.537791
1993 8.676514
1994 9.525407
1995 8.192424
1996 8.622041
1997 8.670823
1998 8.768147
1999 7.907919
2000 8.392275
2001 9.164446
2002 9.813834
2003 8.961799
2004 8.571169
2005 8.096051
2006 8.582421
2007 8.845773
2008 8.787455
2010 8.872551
2011 8.348895
2012 8.502421
2013 8.997153
2014 8.050838
2015 7.657484
2016 8.378217
2017 8.126881
2018 8.079436
2019 8.651915
2020 8.964956
2022 8.545430
ggplot(standard_dev_all_426, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Non-corrected standard deviation of SNOTEL 426 average temperatures for water years 2005-2021

MK & SS for 426 (non-corrected)

sd_mk_426 <- mk.test(standard_dev_all_426$sd_2)
print(sd_mk_426)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_426$sd_2
## z = -1.2749, n = 34, p-value = 0.2023
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  -87.0000000 4550.3333333   -0.1550802
sd_sens_426 <- sens.slope(standard_dev_all_426$sd_2)
print(sd_sens_426)
## 
##  Sen's slope
## 
## data:  standard_dev_all_426$sd_2
## z = -1.2749, n = 34, p-value = 0.2023
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.02882424  0.00714949
## sample estimates:
##  Sen's slope 
## -0.009821212

Corrected

426 Detrended (corrected)

#using the clean culled df:
#average water year temperature
yearly_wy_aver_426_ad <- snotel_426_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_426_ad <- yearly_wy_aver_426_ad %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_426_ad <- daily_wy_aver_426_ad %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp_ad = mean(daily_wy_aver_426_ad$aver_day_temp_ad))
# try to show all years as means. 
daily_wy_aver2_426_ad <-daily_wy_aver_426_ad %>% 
  group_by(waterDay) %>%
  mutate(date_temp_ad = mean(temp_ad))
  
daily_wy_aver2_426_ad$date_temp_ad <- signif(daily_wy_aver2_426_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_426_ad, aes(x = waterDay, y = date_temp_ad))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

426 SS (corrected)

standard_dev_426_ad <- daily_wy_aver_426_ad %>% 
  group_by(waterYear) %>% 
  #filter(waterYear >= 1987 & waterYear <= 2021) %>% 
  mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>% 
  mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_426_ad <- standard_dev_426_ad %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())
standard_dev_all_426_ad <- standard_dev_all_426_ad %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)
standard_dev_all_426_ad %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1987 8.032449
1988 7.469399
1989 8.830085
1990 8.043543
1991 8.571182
1992 7.909178
1993 8.068779
1994 8.814518
1995 7.580772
1996 7.973712
1997 8.044413
1998 8.100725
1999 7.332358
2000 7.718765
2001 8.465805
2002 9.079917
2003 8.257697
2004 7.940383
2005 7.444764
2006 8.582735
2007 8.846403
2008 8.788093
2010 8.872827
2011 8.349023
2012 8.502544
2013 8.997369
2014 8.051074
2015 7.658079
2016 8.378828
2017 8.126888
2018 8.080106
2019 8.652534
2020 8.965133
2022 8.546170
ggplot(standard_dev_all_426_ad, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Morrisey-corrected standard deviation of SNOTEL 426 average temperatures for water years 1986-2021

MK & SS 426 (corrected)

sd_mk_426_ad <- mk.test(standard_dev_all_426_ad$sd_2)
print(sd_mk_426_ad)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_426_ad$sd_2
## z = 1.9272, n = 34, p-value = 0.05396
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  131.0000000 4550.3333333    0.2335116
sd_sens_426_ad <- sens.slope(standard_dev_all_426_ad$sd_2)
print(sd_sens_426_ad)
## 
##  Sen's slope
## 
## data:  standard_dev_all_426_ad$sd_2
## z = 1.9272, n = 34, p-value = 0.05396
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.0008852415  0.0327236747
## sample estimates:
## Sen's slope 
##  0.01613353

Dry Lake 457

Oyler -> Morrisey 7/30/2003

snotel_457 <- SNOTEL_yampa_area %>% 
  filter(site_id == "457")
#str(snotel_457) # check the date, usually a character.  

snotel_457$Date <- as.Date(snotel_457$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.

#THIS WILL CHANGE FOR EACH STATION
snotel_457_clean <- snotel_457 %>% # filter for the timeframe
  filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
  #filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers   
  addWaterYear() %>% 
  mutate(daymonth = format(as.Date(Date), "%d-%m")) %>% 
  na.omit()

#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))

snotel_457_clean <- snotel_457_clean %>% 
  group_by(waterYear)%>% 
  mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers

ggplot(snotel_457_clean, aes(x = Date, y = temperature_mean)) +
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

snotel_457_clean <- snotel_457_clean %>% 
  mutate(temp_diff = abs(temperature_min - temperature_max)) %>% 
  filter(temperature_mean > -40) %>% 
  filter(temperature_mean < 40)
ggplot(snotel_457_clean, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”

# filtering for temperature anomalies
#snotel_457_cull_count <- snotel_457_clean %>% 
#  filter(temperature_min > -40) %>% 
#  count(waterYear)

#snotel_457_cull_count

# filtering for too few observations in a year
snotel_457_cull_count_days <- snotel_457_clean %>% 
  group_by(waterYear) %>% 
  count(waterYear) %>% 
  filter(n < 350)

snotel_457_cull_count_days
## # A tibble: 6 x 2
## # Groups:   waterYear [6]
##   waterYear     n
##       <dbl> <int>
## 1      1985     1
## 2      1994   332
## 3      1996   349
## 4      1997   329
## 5      2003   346
## 6      2021   346
snotel_457_clean_culled <- snotel_457_clean %>% 
  filter(waterYear != "1985" & waterYear != "1986" & waterYear != "1994" & waterYear != "1996" & waterYear != "1997" & waterYear != "2003" & waterYear != "2021")# & waterYear != "1987" & waterYear != "1994" & waterYear != "2002")# & waterYear != "2002" & waterYear != "2016" & waterYear != "2022")# & waterYear != "2017") #%>% 
  #filter(temperature_mean > -49)

ggplot(snotel_457_clean_culled, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

Also 1986

ggplot(snotel_457_clean_culled, aes(x = Date, y = temperature_mean)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

temp_457_xts <- xts(snotel_457_clean_culled$temperature_mean, order.by = snotel_457_clean_culled$Date)

dygraph(temp_457_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 
#snotel_457_clean_culled <- snotel_457_clean_culled %>% 
#  filter(temperature_mean > -30)

#temp_457_xts <- xts(snotel_457_clean_culled$temperature_mean, order.by = snotel_457_clean_culled$Date)

#dygraph(temp_457_xts) %>%
#  dyAxis("y", label = "Daily mean temperature (°C)") 

Dry Lake 457 Oyler -> Morrisey 7/30/2003

snotel_457_adjusted <- snotel_457_clean_culled %>%
  mutate(temp_ad = if_else(Date < "2003-07-30", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))

Non-corrected

457 Detrended

#using the clean culled df:

#average water year temperature

yearly_wy_aver_457 <- snotel_457_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(temperature_mean))

#Average temperature by day for all water years:

daily_wy_aver_457 <- yearly_wy_aver_457 %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

daily_wy_aver_457 <- daily_wy_aver_457 %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(daily_wy_aver_457$aver_day_temp))

# try to show all years as means. 
daily_wy_aver2_457 <-daily_wy_aver_457 %>% 
  group_by(waterDay) %>%
  mutate(date_temp = mean(temperature_mean))
  

daily_wy_aver2_457$date_temp <- signif(daily_wy_aver2_457$date_temp,3) #reduce the sig figs


ggplot(daily_wy_aver2_457, aes(x = waterDay, y = date_temp))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

457 SD

standard_dev_457 <- daily_wy_aver_457 %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

standard_dev_all_457 <- standard_dev_457 %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

standard_dev_all_457 <- standard_dev_all_457 %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

standard_dev_all_457 %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1987 8.892848
1988 9.953110
1989 9.757795
1990 9.176039
1991 9.421986
1992 8.507405
1993 8.878013
1995 8.358076
1998 9.150633
1999 8.440975
2000 8.736877
2001 9.645923
2002 10.040430
2004 8.051569
2005 7.863386
2006 8.768314
2007 9.036800
2008 8.900491
2009 7.976123
2010 8.677359
2011 8.610741
2012 8.805624
2013 9.212073
2014 8.350688
2015 7.879517
2016 8.589405
2017 8.326662
2018 8.323975
2019 8.654312
2020 9.071537
2022 8.555147
ggplot(standard_dev_all_457, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Non-corrected standard deviation of SNOTEL 457 average temperatures for water years 2005-2021

MK & SS for 457 (non-corrected)

sd_mk_457 <- mk.test(standard_dev_all_457$sd_2)
print(sd_mk_457)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_457$sd_2
## z = -2.3795, n = 31, p-value = 0.01734
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## -141.0000000 3461.6666667   -0.3032258
sd_sens_457 <- sens.slope(standard_dev_all_457$sd_2)
print(sd_sens_457)
## 
##  Sen's slope
## 
## data:  standard_dev_all_457$sd_2
## z = -2.3795, n = 31, p-value = 0.01734
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.049787462 -0.005081117
## sample estimates:
## Sen's slope 
## -0.02638274

Corrected

457 Detrended (corrected)

#using the clean culled df:
#average water year temperature
yearly_wy_aver_457_ad <- snotel_457_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_457_ad <- yearly_wy_aver_457_ad %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_457_ad <- daily_wy_aver_457_ad %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp_ad = mean(daily_wy_aver_457_ad$aver_day_temp_ad))
# try to show all years as means. 
daily_wy_aver2_457_ad <-daily_wy_aver_457_ad %>% 
  group_by(waterDay) %>%
  mutate(date_temp_ad = mean(temp_ad))
  
daily_wy_aver2_457_ad$date_temp_ad <- signif(daily_wy_aver2_457_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_457_ad, aes(x = waterDay, y = date_temp_ad))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

457 SS (corrected)

standard_dev_457_ad <- daily_wy_aver_457_ad %>% 
  group_by(waterYear) %>% 
  #filter(waterYear >= 1987 & waterYear <= 2021) %>% 
  mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>% 
  mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_457_ad <- standard_dev_457_ad %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())
standard_dev_all_457_ad <- standard_dev_all_457_ad %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)
standard_dev_all_457_ad %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1987 8.200529
1988 9.184255
1989 9.050691
1990 8.442843
1991 8.739959
1992 7.843690
1993 8.207632
1995 7.712921
1998 8.438042
1999 7.811811
2000 8.022352
2001 8.894507
2002 9.283490
2004 8.051478
2005 7.863308
2006 8.768371
2007 9.036790
2008 8.900641
2009 7.976349
2010 8.677249
2011 8.611094
2012 8.805285
2013 9.211821
2014 8.350803
2015 7.879602
2016 8.589409
2017 8.326313
2018 8.323899
2019 8.654521
2020 9.071250
2022 8.555395
ggplot(standard_dev_all_457_ad, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Morrisey-corrected standard deviation of SNOTEL 457 average temperatures for water years 1986-2021

MK & SS 457 (corrected)

sd_mk_457_ad <- mk.test(standard_dev_all_457_ad$sd_2)
print(sd_mk_457_ad)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_457_ad$sd_2
## z = 0.44191, n = 31, p-value = 0.6586
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 2.700000e+01 3.461667e+03 5.806452e-02
sd_sens_457_ad <- sens.slope(standard_dev_all_457_ad$sd_2)
print(sd_sens_457_ad)
## 
##  Sen's slope
## 
## data:  standard_dev_all_457_ad$sd_2
## z = 0.44191, n = 31, p-value = 0.6586
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01772587  0.02844930
## sample estimates:
## Sen's slope 
## 0.004837867

Elk River 467

Oyler -> Morrisey 8/7/2006

snotel_467 <- SNOTEL_yampa_area %>% 
  filter(site_id == "467")
#str(snotel_467) # check the date, usually a character.  

snotel_467$Date <- as.Date(snotel_467$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.

#THIS WILL CHANGE FOR EACH STATION
snotel_467_clean <- snotel_467 %>% # filter for the timeframe
  filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
  #filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers   
  addWaterYear() %>% 
  mutate(daymonth = format(as.Date(Date), "%d-%m")) %>% 
  na.omit()

#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))

snotel_467_clean <- snotel_467_clean %>% 
  group_by(waterYear)%>% 
  mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers

ggplot(snotel_467_clean, aes(x = Date, y = temperature_mean)) +
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

snotel_467_clean <- snotel_467_clean %>% 
  mutate(temp_diff = abs(temperature_min - temperature_max)) %>% 
  filter(temperature_mean > -40) %>% 
  filter(temperature_mean < 40)
ggplot(snotel_467_clean, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”

# filtering for temperature anomalies
#snotel_467_cull_count <- snotel_467_clean %>% 
#  filter(temperature_min > -40) %>% 
#  count(waterYear)

#snotel_467_cull_count

# filtering for too few observations in a year
snotel_467_cull_count_days <- snotel_467_clean %>% 
  group_by(waterYear) %>% 
  count(waterYear) %>% 
  filter(n < 350)

snotel_467_cull_count_days
## # A tibble: 6 x 2
## # Groups:   waterYear [6]
##   waterYear     n
##       <dbl> <int>
## 1      1985     1
## 2      1994   326
## 3      1998   342
## 4      1999   281
## 5      2001   339
## 6      2021   344
snotel_467_clean_culled <- snotel_467_clean %>% 
  filter(waterYear != "1985" & waterYear != "1986" & waterYear != "1994" & waterYear != "1998" & waterYear != "1999" & waterYear != "2001" & waterYear != "2021")# & waterYear != "2021")# & waterYear != "1987" & waterYear != "1994" & waterYear != "2002")# & waterYear != "2002" & waterYear != "2016" & waterYear != "2022")# & waterYear != "2017") #%>% 
  #filter(temperature_mean > -49)

ggplot(snotel_467_clean_culled, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

Also 1986

ggplot(snotel_467_clean_culled, aes(x = Date, y = temperature_mean)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

temp_467_xts <- xts(snotel_467_clean_culled$temperature_mean, order.by = snotel_467_clean_culled$Date)

dygraph(temp_467_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 
#snotel_467_clean_culled <- snotel_467_clean_culled %>% 
#  filter(temperature_mean > -30)

#temp_467_xts <- xts(snotel_467_clean_culled$temperature_mean, order.by = snotel_467_clean_culled$Date)

#dygraph(temp_467_xts) %>%
#  dyAxis("y", label = "Daily mean temperature (°C)") 

Elk River 467 Oyler -> Morrisey 8/7/2006

snotel_467_adjusted <- snotel_467_clean_culled %>%
  mutate(temp_ad = if_else(Date < "2006-08-07", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))

Non-corrected

467 Detrended

#using the clean culled df:

#average water year temperature

yearly_wy_aver_467 <- snotel_467_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(temperature_mean))

#Average temperature by day for all water years:

daily_wy_aver_467 <- yearly_wy_aver_467 %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

daily_wy_aver_467 <- daily_wy_aver_467 %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(daily_wy_aver_467$aver_day_temp))

# try to show all years as means. 
daily_wy_aver2_467 <-daily_wy_aver_467 %>% 
  group_by(waterDay) %>%
  mutate(date_temp = mean(temperature_mean))
  

daily_wy_aver2_467$date_temp <- signif(daily_wy_aver2_467$date_temp,3) #reduce the sig figs


ggplot(daily_wy_aver2_467, aes(x = waterDay, y = date_temp))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

467 SD

standard_dev_467 <- daily_wy_aver_467 %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

standard_dev_all_467 <- standard_dev_467 %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

standard_dev_all_467 <- standard_dev_all_467 %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

standard_dev_all_467 %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1987 9.103406
1988 10.510726
1989 10.049523
1990 9.547427
1991 9.692602
1992 8.717383
1993 9.114798
1995 8.625842
1996 8.985964
1997 9.119643
2000 8.850761
2002 10.120172
2003 9.161089
2004 8.589463
2005 8.525521
2006 9.564659
2007 9.088368
2008 9.090316
2009 8.102303
2010 8.886760
2011 8.734057
2012 8.902720
2013 9.375328
2014 8.477094
2015 8.044768
2016 8.660871
2017 8.514442
2018 8.465435
2019 8.883368
2020 9.345838
2022 8.902383
ggplot(standard_dev_all_467, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Non-corrected standard deviation of SNOTEL 467 average temperatures for water years 2005-2021

MK & SS for 467 (non-corrected)

sd_mk_467 <- mk.test(standard_dev_all_467$sd_2)
print(sd_mk_467)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_467$sd_2
## z = -2.7194, n = 31, p-value = 0.00654
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## -161.0000000 3461.6666667   -0.3462366
sd_sens_467 <- sens.slope(standard_dev_all_467$sd_2)
print(sd_sens_467)
## 
##  Sen's slope
## 
## data:  standard_dev_all_467$sd_2
## z = -2.7194, n = 31, p-value = 0.00654
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.053188331 -0.007753413
## sample estimates:
## Sen's slope 
## -0.02739628

Corrected

467 Detrended (corrected)

#using the clean culled df:
#average water year temperature
yearly_wy_aver_467_ad <- snotel_467_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_467_ad <- yearly_wy_aver_467_ad %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_467_ad <- daily_wy_aver_467_ad %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp_ad = mean(daily_wy_aver_467_ad$aver_day_temp_ad))
# try to show all years as means. 
daily_wy_aver2_467_ad <-daily_wy_aver_467_ad %>% 
  group_by(waterDay) %>%
  mutate(date_temp_ad = mean(temp_ad))
  
daily_wy_aver2_467_ad$date_temp_ad <- signif(daily_wy_aver2_467_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_467_ad, aes(x = waterDay, y = date_temp_ad))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

467 SS (corrected)

standard_dev_467_ad <- daily_wy_aver_467_ad %>% 
  group_by(waterYear) %>% 
  #filter(waterYear >= 1987 & waterYear <= 2021) %>% 
  mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>% 
  mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_467_ad <- standard_dev_467_ad %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())
standard_dev_all_467_ad <- standard_dev_all_467_ad %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)
standard_dev_all_467_ad %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1987 8.416565
1988 9.721370
1989 9.348058
1990 8.816292
1991 9.016762
1992 8.057564
1993 8.438420
1995 7.963510
1996 8.304891
1997 8.454695
2000 8.126220
2002 9.357756
2003 8.423881
2004 7.952025
2005 7.847462
2006 8.843490
2007 9.088080
2008 9.090597
2009 8.102286
2010 8.886876
2011 8.734381
2012 8.902800
2013 9.375482
2014 8.477382
2015 8.045353
2016 8.661195
2017 8.514535
2018 8.465471
2019 8.883707
2020 9.345830
2022 8.902408
ggplot(standard_dev_all_467_ad, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Morrisey-corrected standard deviation of SNOTEL 467 average temperatures for water years 1986-2021

MK & SS 467 (corrected)

sd_mk_467_ad <- mk.test(standard_dev_all_467_ad$sd_2)
print(sd_mk_467_ad)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_467_ad$sd_2
## z = 0.57788, n = 31, p-value = 0.5633
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 3.500000e+01 3.461667e+03 7.526882e-02
sd_sens_467_ad <- sens.slope(standard_dev_all_467_ad$sd_2)
print(sd_sens_467_ad)
## 
##  Sen's slope
## 
## data:  standard_dev_all_467_ad$sd_2
## z = 0.57788, n = 31, p-value = 0.5633
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01785969  0.02658485
## sample estimates:
## Sen's slope 
## 0.003440178

Lynx Pass 607

Morrisey 5/22/2006

snotel_607 <- SNOTEL_yampa_area %>% 
  filter(site_id == "607")
#str(snotel_607) # check the date, usually a character.  

snotel_607$Date <- as.Date(snotel_607$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.

#THIS WILL CHANGE FOR EACH STATION
snotel_607_clean <- snotel_607 %>% # filter for the timeframe
  filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
  #filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers   
  addWaterYear() %>% 
  mutate(daymonth = format(as.Date(Date), "%d-%m")) %>% 
  na.omit()

#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))

snotel_607_clean <- snotel_607_clean %>% 
  group_by(waterYear)%>% 
  mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers

ggplot(snotel_607_clean, aes(x = Date, y = temperature_mean)) +
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

snotel_607_clean <- snotel_607_clean %>% 
  mutate(temp_diff = abs(temperature_min - temperature_max)) %>% 
  filter(temperature_mean > -40) %>% 
  filter(temperature_mean < 40)
ggplot(snotel_607_clean, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”

# filtering for temperature anomalies
#snotel_607_cull_count <- snotel_607_clean %>% 
#  filter(temperature_min > -40) %>% 
#  count(waterYear)

#snotel_607_cull_count

# filtering for too few observations in a year
snotel_607_cull_count_days <- snotel_607_clean %>% 
  group_by(waterYear) %>% 
  count(waterYear) %>% 
  filter(n < 350)

snotel_607_cull_count_days
## # A tibble: 4 x 2
## # Groups:   waterYear [4]
##   waterYear     n
##       <dbl> <int>
## 1      1985     1
## 2      1994   333
## 3      2019   284
## 4      2020   328
snotel_607_clean_culled <- snotel_607_clean %>% 
  filter(waterYear != "1985" & waterYear != "1986" & waterYear != "1994" & waterYear != "2019" & waterYear != "2020")# & waterYear != "2001" & waterYear != "2021")# & waterYear != "2021")# & waterYear != "1987" & waterYear != "1994" & waterYear != "2002")# & waterYear != "2002" & waterYear != "2016" & waterYear != "2022")# & waterYear != "2017") #%>% 
  #filter(temperature_mean > -49)

ggplot(snotel_607_clean_culled, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

Also 1986

ggplot(snotel_607_clean_culled, aes(x = Date, y = temperature_mean)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

temp_607_xts <- xts(snotel_607_clean_culled$temperature_mean, order.by = snotel_607_clean_culled$Date)

dygraph(temp_607_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 
#snotel_607_clean_culled <- snotel_607_clean_culled %>% 
#  filter(temperature_mean > -30)

#temp_607_xts <- xts(snotel_607_clean_culled$temperature_mean, order.by = snotel_607_clean_culled$Date)

#dygraph(temp_607_xts) %>%
#  dyAxis("y", label = "Daily mean temperature (°C)") 

Lynx Pass 607 Morrisey 5/22/2006

snotel_607_adjusted <- snotel_607_clean_culled %>%
  mutate(temp_ad = if_else(Date < "2006-05-22", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))

Non-corrected

607 Detrended

#using the clean culled df:

#average water year temperature

yearly_wy_aver_607 <- snotel_607_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(temperature_mean))

#Average temperature by day for all water years:

daily_wy_aver_607 <- yearly_wy_aver_607 %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

daily_wy_aver_607 <- daily_wy_aver_607 %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(daily_wy_aver_607$aver_day_temp))

# try to show all years as means. 
daily_wy_aver2_607 <-daily_wy_aver_607 %>% 
  group_by(waterDay) %>%
  mutate(date_temp = mean(temperature_mean))
  

daily_wy_aver2_607$date_temp <- signif(daily_wy_aver2_607$date_temp,3) #reduce the sig figs


ggplot(daily_wy_aver2_607, aes(x = waterDay, y = date_temp))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

607 SD

standard_dev_607 <- daily_wy_aver_607 %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

standard_dev_all_607 <- standard_dev_607 %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

standard_dev_all_607 <- standard_dev_all_607 %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

standard_dev_all_607 %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1987 8.658921
1988 9.691007
1989 9.677093
1990 9.127534
1991 9.267600
1992 8.495617
1993 8.542784
1995 8.427463
1996 8.936524
1997 8.897439
1998 9.097517
1999 8.407095
2000 8.685669
2001 9.559462
2002 10.100422
2003 8.935295
2004 8.581025
2005 8.571462
2006 9.563477
2007 9.027055
2008 8.751128
2009 7.935985
2010 9.022507
2011 8.612122
2012 8.812590
2013 9.410619
2014 8.353556
2015 7.907493
2016 8.522596
2017 8.366193
2018 8.240463
2021 8.972977
2022 8.767136
ggplot(standard_dev_all_607, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Non-corrected standard deviation of SNOTEL 607 average temperatures for water years 2005-2021

MK & SS for 607 (non-corrected)

sd_mk_607 <- mk.test(standard_dev_all_607$sd_2)
print(sd_mk_607)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_607$sd_2
## z = -2.0298, n = 33, p-value = 0.04238
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##        S     varS      tau 
## -132.000 4165.333   -0.250
sd_sens_607 <- sens.slope(standard_dev_all_607$sd_2)
print(sd_sens_607)
## 
##  Sen's slope
## 
## data:  standard_dev_all_607$sd_2
## z = -2.0298, n = 33, p-value = 0.04238
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.036056261 -0.001515865
## sample estimates:
## Sen's slope 
## -0.01788946

Corrected

607 Detrended (corrected)

#using the clean culled df:
#average water year temperature
yearly_wy_aver_607_ad <- snotel_607_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_607_ad <- yearly_wy_aver_607_ad %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_607_ad <- daily_wy_aver_607_ad %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp_ad = mean(daily_wy_aver_607_ad$aver_day_temp_ad))
# try to show all years as means. 
daily_wy_aver2_607_ad <-daily_wy_aver_607_ad %>% 
  group_by(waterDay) %>%
  mutate(date_temp_ad = mean(temp_ad))
  
daily_wy_aver2_607_ad$date_temp_ad <- signif(daily_wy_aver2_607_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_607_ad, aes(x = waterDay, y = date_temp_ad))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

607 SS (corrected)

standard_dev_607_ad <- daily_wy_aver_607_ad %>% 
  group_by(waterYear) %>% 
  #filter(waterYear >= 1987 & waterYear <= 2021) %>% 
  mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>% 
  mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_607_ad <- standard_dev_607_ad %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())
standard_dev_all_607_ad <- standard_dev_all_607_ad %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)
standard_dev_all_607_ad %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1987 8.065805
1988 9.041423
1989 9.118850
1990 8.514101
1991 8.699806
1992 7.927056
1993 7.982533
1995 7.846909
1996 8.320426
1997 8.291259
1998 8.445612
1999 7.829298
2000 8.031514
2001 8.877490
2002 9.409575
2003 8.273822
2004 8.002022
2005 7.966244
2006 8.852307
2007 9.027252
2008 8.750916
2009 7.936108
2010 9.022503
2011 8.612197
2012 8.812993
2013 9.410822
2014 8.353805
2015 7.907862
2016 8.523025
2017 8.366000
2018 8.241006
2021 8.973005
2022 8.767834
ggplot(standard_dev_all_607_ad, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Morrisey-corrected standard deviation of SNOTEL 607 average temperatures for water years 1986-2021

MK & SS 607 (corrected)

sd_mk_607_ad <- mk.test(standard_dev_all_607_ad$sd_2)
print(sd_mk_607_ad)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_607_ad$sd_2
## z = 0.63527, n = 33, p-value = 0.5253
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 4.200000e+01 4.165333e+03 7.954545e-02
sd_sens_607_ad <- sens.slope(standard_dev_all_607_ad$sd_2)
print(sd_sens_607_ad)
## 
##  Sen's slope
## 
## data:  standard_dev_all_607_ad$sd_2
## z = 0.63527, n = 33, p-value = 0.5253
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.008825448  0.026678756
## sample estimates:
## Sen's slope 
## 0.006477733

Rabbit Ears 709

Morrisey 8/7/2006

snotel_709 <- SNOTEL_yampa_area %>% 
  filter(site_id == "709")
#str(snotel_709) # check the date, usually a character.  

snotel_709$Date <- as.Date(snotel_709$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.

#THIS WILL CHANGE FOR EACH STATION
snotel_709_clean <- snotel_709 %>% # filter for the timeframe
  filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
  #filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers   
  addWaterYear() %>% 
  mutate(daymonth = format(as.Date(Date), "%d-%m")) %>% 
  na.omit()

#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))

snotel_709_clean <- snotel_709_clean %>% 
  group_by(waterYear)%>% 
  mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers

ggplot(snotel_709_clean, aes(x = Date, y = temperature_mean)) +
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

snotel_709_clean <- snotel_709_clean %>% 
  mutate(temp_diff = abs(temperature_min - temperature_max)) %>% 
  filter(temperature_mean > -40) %>% 
  filter(temperature_mean < 40)
ggplot(snotel_709_clean, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”

# filtering for temperature anomalies
#snotel_709_cull_count <- snotel_709_clean %>% 
#  filter(temperature_min > -40) %>% 
#  count(waterYear)

#snotel_709_cull_count

# filtering for too few observations in a year
snotel_709_cull_count_days <- snotel_709_clean %>% 
  group_by(waterYear) %>% 
  count(waterYear) %>% 
  filter(n < 350)

snotel_709_cull_count_days
## # A tibble: 14 x 2
## # Groups:   waterYear [14]
##    waterYear     n
##        <dbl> <int>
##  1      1992   318
##  2      1993   293
##  3      1994   332
##  4      1995   169
##  5      2006   348
##  6      2007   311
##  7      2011   292
##  8      2012   308
##  9      2013   300
## 10      2014   283
## 11      2015    44
## 12      2016    18
## 13      2017    34
## 14      2018   147
snotel_709_clean_culled <- snotel_709_clean %>% 
  filter(waterYear > "1995" & waterYear != "2006" & waterYear != "2007" & waterYear != "2011" & waterYear != "2012" & waterYear != "2013" & waterYear != "2014" & waterYear != "2015" & waterYear != "2016" & waterYear != "2017" & waterYear != "2018")# & waterYear != "2002")# & waterYear != "2002" & waterYear != "2016" & waterYear != "2022")# & waterYear != "2017") #%>% 
  #filter(temperature_mean > -49)

ggplot(snotel_709_clean_culled, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

ggplot(snotel_709_clean_culled, aes(x = Date, y = temperature_mean)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

temp_709_xts <- xts(snotel_709_clean_culled$temperature_mean, order.by = snotel_709_clean_culled$Date)

dygraph(temp_709_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 
#snotel_709_clean_culled <- snotel_709_clean_culled %>% 
#  filter(temperature_mean > -30)

#temp_709_xts <- xts(snotel_709_clean_culled$temperature_mean, order.by = snotel_709_clean_culled$Date)

#dygraph(temp_709_xts) %>%
#  dyAxis("y", label = "Daily mean temperature (°C)") 

Rabbit Ears 709 Morrisey 8/7/2006

snotel_709_adjusted <- snotel_709_clean_culled %>%
  mutate(temp_ad = if_else(Date < "2006-08-07", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))

Non-corrected

709 Detrended

#using the clean culled df:

#average water year temperature

yearly_wy_aver_709 <- snotel_709_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(temperature_mean))

#Average temperature by day for all water years:

daily_wy_aver_709 <- yearly_wy_aver_709 %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

daily_wy_aver_709 <- daily_wy_aver_709 %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(daily_wy_aver_709$aver_day_temp))

# try to show all years as means. 
daily_wy_aver2_709 <-daily_wy_aver_709 %>% 
  group_by(waterDay) %>%
  mutate(date_temp = mean(temperature_mean))
  

daily_wy_aver2_709$date_temp <- signif(daily_wy_aver2_709$date_temp,3) #reduce the sig figs


ggplot(daily_wy_aver2_709, aes(x = waterDay, y = date_temp))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

709 SD

standard_dev_709 <- daily_wy_aver_709 %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

standard_dev_all_709 <- standard_dev_709 %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

standard_dev_all_709 <- standard_dev_all_709 %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

standard_dev_all_709 %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1996 8.845436
1997 8.864608
1998 9.079136
1999 8.309577
2000 8.748503
2001 9.426260
2002 9.843776
2003 9.011443
2004 8.578795
2005 8.457080
2008 9.038776
2009 8.123905
2010 8.958458
2019 8.740431
2020 9.196584
2021 9.081772
2022 8.759943
ggplot(standard_dev_all_709, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Non-corrected standard deviation of SNOTEL 709 average temperatures for water years 2005-2021

MK & SS for 709 (non-corrected)

sd_mk_709 <- mk.test(standard_dev_all_709$sd_2)
print(sd_mk_709)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_709$sd_2
## z = 0.041193, n = 17, p-value = 0.9671
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##   2.00000000 589.33333333   0.01470588
sd_sens_709 <- sens.slope(standard_dev_all_709$sd_2)
print(sd_sens_709)
## 
##  Sen's slope
## 
## data:  standard_dev_all_709$sd_2
## z = 0.041193, n = 17, p-value = 0.9671
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.04647212  0.03945218
## sample estimates:
##  Sen's slope 
## 0.0005780099

Corrected

709 Detrended (corrected)

#using the clean culled df:
#average water year temperature
yearly_wy_aver_709_ad <- snotel_709_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_709_ad <- yearly_wy_aver_709_ad %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_709_ad <- daily_wy_aver_709_ad %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp_ad = mean(daily_wy_aver_709_ad$aver_day_temp_ad))
# try to show all years as means. 
daily_wy_aver2_709_ad <-daily_wy_aver_709_ad %>% 
  group_by(waterDay) %>%
  mutate(date_temp_ad = mean(temp_ad))
  
daily_wy_aver2_709_ad$date_temp_ad <- signif(daily_wy_aver2_709_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_709_ad, aes(x = waterDay, y = date_temp_ad))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

709 SS (corrected)

standard_dev_709_ad <- daily_wy_aver_709_ad %>% 
  group_by(waterYear) %>% 
  #filter(waterYear >= 1987 & waterYear <= 2021) %>% 
  mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>% 
  mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_709_ad <- standard_dev_709_ad %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())
standard_dev_all_709_ad <- standard_dev_all_709_ad %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)
standard_dev_all_709_ad %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1996 8.213076
1997 8.255603
1998 8.412792
1999 7.726856
2000 8.073670
2001 8.727586
2002 9.142221
2003 8.322566
2004 7.979937
2005 7.819602
2008 9.038980
2009 8.124255
2010 8.958543
2019 8.740727
2020 9.196244
2021 9.081146
2022 8.760212
ggplot(standard_dev_all_709_ad, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Morrisey-corrected standard deviation of SNOTEL 709 average temperatures for water years 1986-2021

MK & SS 709 (corrected)

sd_mk_709_ad <- mk.test(standard_dev_all_709_ad$sd_2)
print(sd_mk_709_ad)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_709_ad$sd_2
## z = 1.9361, n = 17, p-value = 0.05286
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##           S        varS         tau 
##  48.0000000 589.3333333   0.3529412
sd_sens_709_ad <- sens.slope(standard_dev_all_709_ad$sd_2)
print(sd_sens_709_ad)
## 
##  Sen's slope
## 
## data:  standard_dev_all_709_ad$sd_2
## z = 1.9361, n = 17, p-value = 0.05286
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.006786111  0.102902003
## sample estimates:
## Sen's slope 
##  0.05101403

Ripple Creek 717

Oyler -> Morrisey 8/7/2006

snotel_717 <- SNOTEL_yampa_area %>% 
  filter(site_id == "717")
#str(snotel_717) # check the date, usually a character.  

snotel_717$Date <- as.Date(snotel_717$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.

#THIS WILL CHANGE FOR EACH STATION
snotel_717_clean <- snotel_717 %>% # filter for the timeframe
  filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
  #filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers   
  addWaterYear() %>% 
  mutate(daymonth = format(as.Date(Date), "%d-%m")) %>% 
  na.omit()

#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))

snotel_717_clean <- snotel_717_clean %>% 
  group_by(waterYear)%>% 
  mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers

ggplot(snotel_717_clean, aes(x = Date, y = temperature_mean)) +
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

snotel_717_clean <- snotel_717_clean %>% 
  mutate(temp_diff = abs(temperature_min - temperature_max)) %>% 
  filter(temperature_mean > -40) %>% 
  filter(temperature_mean < 40)
ggplot(snotel_717_clean, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”

# filtering for temperature anomalies
#snotel_717_cull_count <- snotel_717_clean %>% 
#  filter(temperature_min > -40) %>% 
#  count(waterYear)

#snotel_717_cull_count

# filtering for too few observations in a year
snotel_717_cull_count_days <- snotel_717_clean %>% 
  group_by(waterYear) %>% 
  count(waterYear) %>% 
  filter(n < 350)

snotel_717_cull_count_days
## # A tibble: 2 x 2
## # Groups:   waterYear [2]
##   waterYear     n
##       <dbl> <int>
## 1      1992   311
## 2      1993   330
snotel_717_clean_culled <- snotel_717_clean %>% 
  filter(waterYear != "1992" & waterYear != "1993")# & waterYear != "2007" & waterYear != "2011" & waterYear != "2012" & waterYear != "2013" & waterYear != "2014" & waterYear != "2015" & waterYear != "2016" & waterYear != "2017" & waterYear != "2018")# & waterYear != "2002")# & waterYear != "2002" & waterYear != "2016" & waterYear != "2022")# & waterYear != "2017") #%>% 
  #filter(temperature_mean > -49)

ggplot(snotel_717_clean_culled, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

ggplot(snotel_717_clean_culled, aes(x = Date, y = temperature_mean)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

temp_717_xts <- xts(snotel_717_clean_culled$temperature_mean, order.by = snotel_717_clean_culled$Date)

dygraph(temp_717_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 
#snotel_717_clean_culled <- snotel_717_clean_culled %>% 
#  filter(temperature_mean > -30)

#temp_717_xts <- xts(snotel_717_clean_culled$temperature_mean, order.by = snotel_717_clean_culled$Date)

#dygraph(temp_717_xts) %>%
#  dyAxis("y", label = "Daily mean temperature (°C)") 

Ripple Creek 717 Oyler -> Morrisey 8/7/2006

snotel_717_adjusted <- snotel_717_clean_culled %>%
  mutate(temp_ad = if_else(Date < "2006-08-07", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))

Non-corrected

717 Detrended

#using the clean culled df:

#average water year temperature

yearly_wy_aver_717 <- snotel_717_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(temperature_mean))

#Average temperature by day for all water years:

daily_wy_aver_717 <- yearly_wy_aver_717 %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

daily_wy_aver_717 <- daily_wy_aver_717 %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(daily_wy_aver_717$aver_day_temp))

# try to show all years as means. 
daily_wy_aver2_717 <-daily_wy_aver_717 %>% 
  group_by(waterDay) %>%
  mutate(date_temp = mean(temperature_mean))
  

daily_wy_aver2_717$date_temp <- signif(daily_wy_aver2_717$date_temp,3) #reduce the sig figs


ggplot(daily_wy_aver2_717, aes(x = waterDay, y = date_temp))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

717 SD

standard_dev_717 <- daily_wy_aver_717 %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

standard_dev_all_717 <- standard_dev_717 %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

standard_dev_all_717 <- standard_dev_all_717 %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

standard_dev_all_717 %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1987 8.803468
1988 9.707652
1989 9.270222
1990 9.163954
1991 9.043756
1994 9.475534
1995 8.216129
1996 8.766674
1997 8.761084
1998 8.941067
1999 8.106802
2000 8.669178
2001 9.209650
2002 9.676307
2003 9.008971
2004 8.362946
2005 8.373378
2006 9.248617
2007 8.800246
2008 8.921082
2009 7.998549
2010 8.681956
2011 8.605749
2012 8.582569
2013 9.046302
2014 8.329385
2015 7.737981
2016 8.402920
2017 8.250374
2018 8.202730
2019 8.544093
2020 8.992894
2021 8.957109
2022 8.519914
ggplot(standard_dev_all_717, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Non-corrected standard deviation of SNOTEL 717 average temperatures for water years 2005-2021

MK & SS for 717 (non-corrected)

sd_mk_717 <- mk.test(standard_dev_all_717$sd_2)
print(sd_mk_717)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_717$sd_2
## z = -2.6684, n = 34, p-value = 0.007621
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## -181.0000000 4550.3333333   -0.3226381
sd_sens_717 <- sens.slope(standard_dev_all_717$sd_2)
print(sd_sens_717)
## 
##  Sen's slope
## 
## data:  standard_dev_all_717$sd_2
## z = -2.6684, n = 34, p-value = 0.007621
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.039224947 -0.006265457
## sample estimates:
## Sen's slope 
## -0.02259117

Corrected

717 Detrended (corrected)

#using the clean culled df:
#average water year temperature
yearly_wy_aver_717_ad <- snotel_717_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_717_ad <- yearly_wy_aver_717_ad %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_717_ad <- daily_wy_aver_717_ad %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp_ad = mean(daily_wy_aver_717_ad$aver_day_temp_ad))
# try to show all years as means. 
daily_wy_aver2_717_ad <-daily_wy_aver_717_ad %>% 
  group_by(waterDay) %>%
  mutate(date_temp_ad = mean(temp_ad))
  
daily_wy_aver2_717_ad$date_temp_ad <- signif(daily_wy_aver2_717_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_717_ad, aes(x = waterDay, y = date_temp_ad))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

717 SS (corrected)

standard_dev_717_ad <- daily_wy_aver_717_ad %>% 
  group_by(waterYear) %>% 
  #filter(waterYear >= 1987 & waterYear <= 2021) %>% 
  mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>% 
  mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_717_ad <- standard_dev_717_ad %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())
standard_dev_all_717_ad <- standard_dev_all_717_ad %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)
standard_dev_all_717_ad %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1987 8.203857
1988 9.063074
1989 8.683957
1990 8.529404
1991 8.477488
1994 8.832494
1995 7.643847
1996 8.162790
1997 8.181297
1998 8.315921
1999 7.566445
2000 8.038773
2001 8.563752
2002 9.013927
2003 8.355528
2004 7.809021
2005 7.779504
2006 8.600879
2007 8.799659
2008 8.920535
2009 7.997923
2010 8.681418
2011 8.605311
2012 8.581992
2013 9.045623
2014 8.328989
2015 7.737670
2016 8.402569
2017 8.249042
2018 8.202635
2019 8.543763
2020 8.992050
2021 8.956530
2022 8.520306
ggplot(standard_dev_all_717_ad, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Morrisey-corrected standard deviation of SNOTEL 717 average temperatures for water years 1986-2021

MK & SS 717 (corrected)

sd_mk_717_ad <- mk.test(standard_dev_all_717_ad$sd_2)
print(sd_mk_717_ad)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_717_ad$sd_2
## z = 0.56333, n = 34, p-value = 0.5732
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 3.900000e+01 4.550333e+03 6.951872e-02
sd_sens_717_ad <- sens.slope(standard_dev_all_717_ad$sd_2)
print(sd_sens_717_ad)
## 
##  Sen's slope
## 
## data:  standard_dev_all_717_ad$sd_2
## z = 0.56333, n = 34, p-value = 0.5732
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01114957  0.02335425
## sample estimates:
## Sen's slope 
## 0.004691148

Tower 825

Original 8/18/2004

snotel_825 <- SNOTEL_yampa_area %>% 
  filter(site_id == "825")
#str(snotel_825) # check the date, usually a character.  

snotel_825$Date <- as.Date(snotel_825$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.

#THIS WILL CHANGE FOR EACH STATION
snotel_825_clean <- snotel_825 %>% # filter for the timeframe
  filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
  #filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers   
  addWaterYear() %>% 
  mutate(daymonth = format(as.Date(Date), "%d-%m")) %>% 
  na.omit()

#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))

snotel_825_clean <- snotel_825_clean %>% 
  group_by(waterYear)%>% 
  mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers

ggplot(snotel_825_clean, aes(x = Date, y = temperature_mean)) +
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

snotel_825_clean <- snotel_825_clean %>% 
  mutate(temp_diff = abs(temperature_min - temperature_max)) %>% 
  filter(temperature_mean > -40) %>% 
  filter(temperature_mean < 30)
ggplot(snotel_825_clean, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”

# filtering for temperature anomalies
#snotel_825_cull_count <- snotel_825_clean %>% 
#  filter(temperature_min > -40) %>% 
#  count(waterYear)

#snotel_825_cull_count

# filtering for too few observations in a year
snotel_825_cull_count_days <- snotel_825_clean %>% 
  group_by(waterYear) %>% 
  count(waterYear) %>% 
  filter(n < 350)

snotel_825_cull_count_days
## # A tibble: 9 x 2
## # Groups:   waterYear [9]
##   waterYear     n
##       <dbl> <int>
## 1      1985     1
## 2      1987   348
## 3      1994   301
## 4      1995   347
## 5      1997   322
## 6      1998   343
## 7      1999   348
## 8      2002   327
## 9      2008   325
snotel_825_clean_culled <- snotel_825_clean %>% 
  filter(waterYear > "1987" & waterYear != "1994" & waterYear != "1995" & waterYear != "1997" & waterYear != "1998" & waterYear != "1999" & waterYear != "2002" & waterYear != "2008")# & waterYear != "2017" & waterYear != "2018")# & waterYear != "2002")# & waterYear != "2002" & waterYear != "2016" & waterYear != "2022")# & waterYear != "2017") #%>% 
  #filter(temperature_mean > -49)

ggplot(snotel_825_clean_culled, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

ggplot(snotel_825_clean_culled, aes(x = Date, y = temperature_mean)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

temp_825_xts <- xts(snotel_825_clean_culled$temperature_mean, order.by = snotel_825_clean_culled$Date)

dygraph(temp_825_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 
#snotel_825_clean_culled <- snotel_825_clean_culled %>% 
#  filter(temperature_mean > -30)

#temp_825_xts <- xts(snotel_825_clean_culled$temperature_mean, order.by = snotel_825_clean_culled$Date)

#dygraph(temp_825_xts) %>%
#  dyAxis("y", label = "Daily mean temperature (°C)") 

Tower 825 Original 8/18/2004

snotel_825_adjusted <- snotel_825_clean_culled %>%
  mutate(temp_ad = if_else(Date < "2004-08-18", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))

Non-corrected

825 Detrended

#using the clean culled df:

#average water year temperature

yearly_wy_aver_825 <- snotel_825_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(temperature_mean))

#Average temperature by day for all water years:

daily_wy_aver_825 <- yearly_wy_aver_825 %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

daily_wy_aver_825 <- daily_wy_aver_825 %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(daily_wy_aver_825$aver_day_temp))

# try to show all years as means. 
daily_wy_aver2_825 <-daily_wy_aver_825 %>% 
  group_by(waterDay) %>%
  mutate(date_temp = mean(temperature_mean))
  

daily_wy_aver2_825$date_temp <- signif(daily_wy_aver2_825$date_temp,3) #reduce the sig figs


ggplot(daily_wy_aver2_825, aes(x = waterDay, y = date_temp))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

825 SD

standard_dev_825 <- daily_wy_aver_825 %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

standard_dev_all_825 <- standard_dev_825 %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

standard_dev_all_825 <- standard_dev_all_825 %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

standard_dev_all_825 %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1988 10.158111
1989 9.521283
1990 9.387663
1991 9.495098
1992 8.468342
1993 8.686860
1996 9.236157
2000 8.951568
2001 9.626977
2003 9.317887
2004 8.794074
2005 8.168991
2006 9.046952
2007 9.112464
2009 8.322586
2010 8.961897
2011 8.871929
2012 8.997194
2013 9.382743
2014 8.649080
2015 8.104451
2016 8.792359
2017 8.593127
2018 8.634239
2019 8.862025
2020 9.328094
2021 9.270956
2022 9.124264
ggplot(standard_dev_all_825, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Non-corrected standard deviation of SNOTEL 825 average temperatures for water years 2005-2021

MK & SS for 825 (non-corrected)

sd_mk_825 <- mk.test(standard_dev_all_825$sd_2)
print(sd_mk_825)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_825$sd_2
## z = -1.7188, n = 28, p-value = 0.08565
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  -88.0000000 2562.0000000   -0.2328042
sd_sens_825 <- sens.slope(standard_dev_all_825$sd_2)
print(sd_sens_825)
## 
##  Sen's slope
## 
## data:  standard_dev_all_825$sd_2
## z = -1.7188, n = 28, p-value = 0.08565
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.04047480  0.00485368
## sample estimates:
## Sen's slope 
## -0.02125455

Corrected

825 Detrended (corrected)

#using the clean culled df:
#average water year temperature
yearly_wy_aver_825_ad <- snotel_825_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_825_ad <- yearly_wy_aver_825_ad %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_825_ad <- daily_wy_aver_825_ad %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp_ad = mean(daily_wy_aver_825_ad$aver_day_temp_ad))
# try to show all years as means. 
daily_wy_aver2_825_ad <-daily_wy_aver_825_ad %>% 
  group_by(waterDay) %>%
  mutate(date_temp_ad = mean(temp_ad))
  
daily_wy_aver2_825_ad$date_temp_ad <- signif(daily_wy_aver2_825_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_825_ad, aes(x = waterDay, y = date_temp_ad))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

825 SS (corrected)

standard_dev_825_ad <- daily_wy_aver_825_ad %>% 
  group_by(waterYear) %>% 
  #filter(waterYear >= 1987 & waterYear <= 2021) %>% 
  mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>% 
  mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_825_ad <- standard_dev_825_ad %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())
standard_dev_all_825_ad <- standard_dev_all_825_ad %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)
standard_dev_all_825_ad %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1988 9.527196
1989 8.967107
1990 8.793798
1991 8.946793
1992 7.952054
1993 8.156645
1996 8.650702
2000 8.338506
2001 8.985410
2003 8.675745
2004 8.209448
2005 8.172186
2006 9.049660
2007 9.115791
2009 8.324666
2010 8.964371
2011 8.875058
2012 9.001462
2013 9.386059
2014 8.652436
2015 8.107386
2016 8.796286
2017 8.596142
2018 8.637340
2019 8.865259
2020 9.331765
2021 9.273785
2022 9.127176
ggplot(standard_dev_all_825_ad, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Morrisey-corrected standard deviation of SNOTEL 825 average temperatures for water years 1986-2021

MK & SS 825 (corrected)

sd_mk_825_ad <- mk.test(standard_dev_all_825_ad$sd_2)
print(sd_mk_825_ad)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_825_ad$sd_2
## z = 1.0076, n = 28, p-value = 0.3137
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##   52.0000000 2562.0000000    0.1375661
sd_sens_825_ad <- sens.slope(standard_dev_all_825_ad$sd_2)
print(sd_sens_825_ad)
## 
##  Sen's slope
## 
## data:  standard_dev_all_825_ad$sd_2
## z = 1.0076, n = 28, p-value = 0.3137
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.009639429  0.034852064
## sample estimates:
## Sen's slope 
##  0.01232876

Trapper Lake 827 Original 12/13/2004

Trapper Lake 827

Original 12/13/2004

snotel_827 <- SNOTEL_yampa_area %>% 
  filter(site_id == "827")
#str(snotel_827) # check the date, usually a character.  

snotel_827$Date <- as.Date(snotel_827$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.

#THIS WILL CHANGE FOR EACH STATION
snotel_827_clean <- snotel_827 %>% # filter for the timeframe
  filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
  #filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers   
  addWaterYear() %>% 
  mutate(daymonth = format(as.Date(Date), "%d-%m")) %>% 
  na.omit()

#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))

snotel_827_clean <- snotel_827_clean %>% 
  group_by(waterYear)%>% 
  mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers

ggplot(snotel_827_clean, aes(x = Date, y = temperature_mean)) +
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

snotel_827_clean <- snotel_827_clean %>% 
  mutate(temp_diff = abs(temperature_min - temperature_max)) %>% 
  filter(temperature_mean > -40) %>% 
  filter(temperature_mean < 30)
ggplot(snotel_827_clean, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”

# filtering for temperature anomalies
#snotel_827_cull_count <- snotel_827_clean %>% 
#  filter(temperature_min > -40) %>% 
#  count(waterYear)

#snotel_827_cull_count

# filtering for too few observations in a year
snotel_827_cull_count_days <- snotel_827_clean %>% 
  group_by(waterYear) %>% 
  count(waterYear) %>% 
  filter(n < 350)

snotel_827_cull_count_days
## # A tibble: 18 x 2
## # Groups:   waterYear [18]
##    waterYear     n
##        <dbl> <int>
##  1      1986   331
##  2      1994   179
##  3      1995   323
##  4      1996    53
##  5      1997   278
##  6      1999   341
##  7      2001   348
##  8      2002   310
##  9      2003   340
## 10      2004   342
## 11      2009   342
## 12      2010   346
## 13      2011   326
## 14      2013   236
## 15      2014   241
## 16      2016   321
## 17      2017   219
## 18      2019   332
snotel_827_clean_culled <- snotel_827_clean %>% 
  filter(waterYear != "1986" & waterYear != "1994" & waterYear != "1995" & waterYear != "1996" & waterYear != "1997" & waterYear != "1999" & waterYear != "2001" & waterYear != "2002" & waterYear != "2003" & waterYear != "2004" & waterYear != "2009" & waterYear != "2010" & waterYear != "2011" & waterYear != "2013" & waterYear != "2014" & waterYear != "2016" & waterYear != "2017" & waterYear != "2019") #%>% 
  #filter(temperature_mean > -49)

ggplot(snotel_827_clean_culled, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

ggplot(snotel_827_clean_culled, aes(x = Date, y = temperature_mean)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

temp_827_xts <- xts(snotel_827_clean_culled$temperature_mean, order.by = snotel_827_clean_culled$Date)

dygraph(temp_827_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 
#snotel_827_clean_culled <- snotel_827_clean_culled %>% 
#  filter(temperature_mean > -30)

#temp_827_xts <- xts(snotel_827_clean_culled$temperature_mean, order.by = snotel_827_clean_culled$Date)

#dygraph(temp_827_xts) %>%
#  dyAxis("y", label = "Daily mean temperature (°C)") 

Trapper Lake 827

Original 12/13/2004

snotel_827_adjusted <- snotel_827_clean_culled %>%
  mutate(temp_ad = if_else(Date < "2004-12-13", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))

Non-corrected

827 Detrended

#using the clean culled df:

#average water year temperature

yearly_wy_aver_827 <- snotel_827_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(temperature_mean))

#Average temperature by day for all water years:

daily_wy_aver_827 <- yearly_wy_aver_827 %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

daily_wy_aver_827 <- daily_wy_aver_827 %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(daily_wy_aver_827$aver_day_temp))

# try to show all years as means. 
daily_wy_aver2_827 <-daily_wy_aver_827 %>% 
  group_by(waterDay) %>%
  mutate(date_temp = mean(temperature_mean))
  

daily_wy_aver2_827$date_temp <- signif(daily_wy_aver2_827$date_temp,3) #reduce the sig figs


ggplot(daily_wy_aver2_827, aes(x = waterDay, y = date_temp))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

827 SD

standard_dev_827 <- daily_wy_aver_827 %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

standard_dev_all_827 <- standard_dev_827 %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

standard_dev_all_827 <- standard_dev_all_827 %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

standard_dev_all_827 %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1987 8.785352
1988 9.643636
1989 9.399426
1990 9.017517
1991 9.053817
1992 8.461780
1993 8.592203
1998 9.226580
2000 8.675610
2005 7.901837
2006 8.811208
2007 8.956871
2008 9.053611
2012 8.641372
2015 7.855779
2018 8.196461
2020 8.933495
2021 8.790312
2022 8.662496
ggplot(standard_dev_all_827, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Non-corrected standard deviation of SNOTEL 827 average temperatures for water years 2005-2021

MK & SS for 827 (non-corrected)

sd_mk_827 <- mk.test(standard_dev_all_827$sd_2)
print(sd_mk_827)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_827$sd_2
## z = -1.7493, n = 19, p-value = 0.08024
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##           S        varS         tau 
## -51.0000000 817.0000000  -0.2982456
sd_sens_827 <- sens.slope(standard_dev_all_827$sd_2)
print(sd_sens_827)
## 
##  Sen's slope
## 
## data:  standard_dev_all_827$sd_2
## z = -1.7493, n = 19, p-value = 0.08024
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.06891402  0.00422489
## sample estimates:
## Sen's slope 
## -0.03456914

Corrected

827 Detrended (corrected)

#using the clean culled df:
#average water year temperature
yearly_wy_aver_827_ad <- snotel_827_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_827_ad <- yearly_wy_aver_827_ad %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_827_ad <- daily_wy_aver_827_ad %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp_ad = mean(daily_wy_aver_827_ad$aver_day_temp_ad))
# try to show all years as means. 
daily_wy_aver2_827_ad <-daily_wy_aver_827_ad %>% 
  group_by(waterDay) %>%
  mutate(date_temp_ad = mean(temp_ad))
  
daily_wy_aver2_827_ad$date_temp_ad <- signif(daily_wy_aver2_827_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_827_ad, aes(x = waterDay, y = date_temp_ad))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

827 SS (corrected)

standard_dev_827_ad <- daily_wy_aver_827_ad %>% 
  group_by(waterYear) %>% 
  #filter(waterYear >= 1987 & waterYear <= 2021) %>% 
  mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>% 
  mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_827_ad <- standard_dev_827_ad %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())
standard_dev_all_827_ad <- standard_dev_all_827_ad %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)
standard_dev_all_827_ad %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1987 8.205336
1988 9.015708
1989 8.826659
1990 8.399945
1991 8.491804
1992 7.896678
1993 8.029814
1998 8.583994
2000 8.036502
2005 7.717045
2006 8.811366
2007 8.957047
2008 9.053812
2012 8.641373
2015 7.856092
2018 8.196931
2020 8.933512
2021 8.790758
2022 8.663052
ggplot(standard_dev_all_827_ad, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Morrisey-corrected standard deviation of SNOTEL 827 average temperatures for water years 1986-2021

MK & SS 827 (corrected)

sd_mk_827_ad <- mk.test(standard_dev_all_827_ad$sd_2)
print(sd_mk_827_ad)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_827_ad$sd_2
## z = 0.41983, n = 19, p-value = 0.6746
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  13.00000000 817.00000000   0.07602339
sd_sens_827_ad <- sens.slope(standard_dev_all_827_ad$sd_2)
print(sd_sens_827_ad)
## 
##  Sen's slope
## 
## data:  standard_dev_all_827_ad$sd_2
## z = 0.41983, n = 19, p-value = 0.6746
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.03007505  0.05895184
## sample estimates:
## Sen's slope 
##    0.012232

Note: figure captions are incorrect.